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Abstract 



Jj^ We describe how to consistently incorporate solar model uncertainties, along with experimental 
^ errors and correlations, when analyzing solar neutrino data to derive confidence limits on parameter 



space for proposed solutions of the solar neutrino problem. Our work resolves ambiguities and 
inconsistencies in the previous literature. As an application of our methods we calculate the masses 



X 

- - and mixing angles allowed by the current data for the proposed MSW solution using both Bayesian 
and frequentist methods, allowing purely for solar model flux variations, to compare with previous 
work. We consider the effects of including metal diffusion in the solar models and also discuss 
implications for future experiments. 
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1. Introduction 

As more experimental information has become available and theorists have converged on "new- 
physics" explanations |IJ, |3|, |], |5], ^ [7| of the solar neutrino problem there has been interest 
in incorporating the error budget of solar models into analyses of the data. This has proceeded in 
stages. The first incorporation of solar model flux uncertainties did not take experimental correlations 
into account |TJ. A subsequent analysis rectified this problem, but then did not properly account for 
neutrino flux correlations f|. Most recently, a detailed analysis has been performed which has largely 
resolved this latter problem by an improved approximation for solar model uncertainties, a correct 
accounting for experimental correlations, as well as a careful examination of such effects as MSW |§ 
mixing in the Earth in order to derive allowed regions of mass and mixing angle [10]. Nevertheless, 
the general applicability of the approximations used there to model solar model uncertainties is not 
guaranteed. In addition, the determination of confidence limits and allowed regions of parameter 
space uses a non-standard statistical analysis. 

Now that it appears that the gallium results are stable and that no new significant experimental 
light is likely to be shed on the problem until the gallium experiments have been checked with 
neutrino sources (GALLEX is scheduled for 'calibration' in June 1994) or the next generation of 
detectors comes on line in 3-4 years, there is time to consider a comprehensive, consistent statistical 
analysis, vis a vis neutrino-based solutions of the solar neutrino problem. (Neutrino, rather than solar 
model, based solutions are now strongly indicated by the present data, even without including the 
Homestake results! ^|, 0]) Such an analysis is the purpose of the present paper. We shall demonstrate 
a technique which treats known solar model uncertainties in a computationally simple fashion, and 
then describe how to incorporate the existing experimental information in order to derive confidence 
limits on neutrino masses and mixing angles which have a well-defined statistical meaning. In the 
approximation in which all solar model uncertainties can be parameterized in terms of the neutrino 
flux uncertainties, this technique yields allowed regions in parameter space which can be compared 
with previous results. 

The determination of allowed regions requires four distinct parts: (1) a calculation of solar model 
uncertainties, (2) a model of neutrino transport and detection probabilities, (3) a determination of 
experimental uncertainties and correlations, and finally (4) a well defined statistical procedure for 
comparing predictions and observations. 

The outline of the paper is as follows. We first describe solar model uncertainties gleaned from 
Monte Carlo studies of solar models. We demonstrate that the essential information about this 
type of solar model uncertainty is contained in the neutrino flux correlation matrix, which can be 
calculated either directly using the solar models themselves, or else using a simple but well defined 
approximation. Next we demonstrate how to translate these flux correlations into an experimental 
covariance matrix necessary to properly incorporate the experimental error budget. Following this 
we describe, for both MSW and vacuum oscillations, how one derives survival probabilities following 
transport through the sun and earth. Finally, we describe how to consistently derive allowed regions 
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for neutrino parameter space using well defined statistical probes in a way which avoids problems 
with past analyses, and discuss the meaning of our results for future experiments and particle physics 
models. 

It is worth emphasizing in advance that by outlining a well defined statistical procedure for 
comparing theory and observation we do not necessarily subscribe to the view that only statistical 
solar model uncertainties are relevant, or even that they may be the most important uncertainties. 
It is quite possible, indeed perhaps likely, that systematic solar model uncertainties, due primarily 
to the introduction of new physics into the model (such as heavy metal diffusion — see below) could 
shift the entire allowed range of model parameters determined by the methods we describe here. 
Nevertheless, as the standard solar model gets more complete, this will be less likely. Our purpose 
here is to define a consistent and correct procedure which may be applied as both the data and the 
theoretical models improve. 

2. Solar Model Uncertainties 

Comprehensive estimates of the present solar model uncertainties have been made by Bahcall and 
collaborators fTl , E| , who performed detailed Monte Carlo analyses of the neutrino fluxes that result 



when solar model input parameters are varied over their allowed ranges. Since calculating many full 
solar models can prove cumbersome in terms of computing time, it is useful to have a reliable and 
efficient approximation scheme which reproduces the results of such a calculation. Several schemes 
have been proposed which account for the variation in the total flux of neutrinos, which is in general 
the major source of uncertainty (from solar physics) in the prediction of the experimental rates. 

Two different approaches have been applied to this problem. The first involves simplifying the 
solar model parameter space, an example of which we will call "the T c approach" ||. Here the fluxes 
are parameterized by a single (solar model output) variable - the core temperature of the sun: T c . 
The temperature dependence of the various fluxes are derived from the scatter plots of flux vs T c 
from solar model Monte Carlo calculations. (See e.g. figures 6.2 and 6.3 in ref ||12|| .) Approximating 
the temperature dependence of the fluxes by power laws in T c specifies the flux distributions, with 
the error in T c determined so as to give the appropriate uncertainties in the fluxes. 

While the scatter plots indicate that the relationship between the neutrino flux and T c can be 
approximately described by a simple power law, this relationship is only approximate and there 
remains a significant width to the straight line that would describe a perfect power law dependence. 
Because of this width, these plots do not indicate how the various fluxes are correlated. For example, 
a solar model with a higher T c may correspond to an increased pp flux, but little or no corresponding 
decrease in the 8 B flux. The T c method, based on only one parameter, of course produces totally 
(anti-) correlated uncertainties for the neutrino fluxes while the solar model flux uncertainties exhibit 
a wide range of correlations. The differences in the correlations for the T c parameterization and the 
full solar model Monte Carlo are a reflection of the scatter in the plots of [|T!J . By overestimating the 



correlations, the T c approach tends to underestimate the size of the allowed parameter region for a 
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given confidence level. Thus while a T c parameterization can be a useful tool in some instances, it is 
not appropriate for calculating solar model uncertainties ||. 

An updated version of this method [nj includes not only T c but also cross section uncertainties in 



the form of two extra parameters, chosen from among the (nuclear cross section) input parameters to 
the solar models. This allows this method to be tuned to more closely approximate the full solar model 
correlations [ ID| . Any method which reproduces the flux correlation matrix can correctly include the 



solar model errors, so this updated method and the "power law method" (described below) should 
agree on the statistical content of the solar model uncertainty. However, the applicability of this 
method, including the determination of which combination works, and which T c uncertainty to use 
can only strictly be determined after the fact by explicitly utilizing the detailed results of the full 
solar model Monte Carlo calculations. 

An alternative approach, proposed earlier |], parameterizes the solar model uncertainties in terms 
of the logarithmic derivatives of the fluxes with respect to the 10 solar model input parameters. It 



was shown in |L2[ that for small variations in the input parameters the neutrino fluxes, 0, can be 
expressed as 

k 

where a^jt is the logarithmic partial derivative of <ftj (j = pp, pep, hep, 7 Be, 8 B, 13 N, 15 and 17 F) with 
respect to the input parameter IV The solar model flux uncertainties can thereafter be obtained from 
a Monte Carlo procedure assuming Gaussian distributions for the input parameters, as described in 
more detail in This method has a firm basis in describing the errors in the output function (solar 
model fluxes) in terms of the errors in the input parameters, and the ctj^ are readily available ||12|| . 
We shall refer to this as "the power law approach" . 

From a Monte Carlo analysis using this approach, we obtain an estimate of the theoretical uncer- 
tainties in the predicted fluxes for each species. This allows us to determine the correlations between 
the various fluxes, which will be important for computing correlations between the rates predicted 
for different detectors. The elements of the covariance matrix for the various fluxes (4>j) are given by 

m 

V jk = U<l>j-fa)(<l>k-4>k)) (2) 



where the angled brackets indicate an average over the solar models and <fi = (0). To display the 
correlations we present the correlation matrix, whose elements are given by 

Pjk = -r~, (3) 

where <jj = yJVjj is the standard deviation in <pj. Note that in the correlation matrix the diagonal 
elements are 1 by definition and off-diagonal elements are equal to (— )1 in the limit of perfect 
(anti-)correlation. 
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It is important to note that the partial derivatives in (|I|) were determined via the solar model 
with fixed solar luminosity. Even though the power law approach does not explicitly enforce such a 
constraint, the use of the relation ([!]) will result in a covariance matrix equal to that from the fully 
self-consistent solar model Monte Carlo calculations. This was implicitly exploited in our previous 
work P], |nj and is explicitly demonstrated in tables |3| and f| which compare the covariance matrices 



for the two approaches. The agreement between our Monte Carlo calculation and the 1,000 models 
of Bahcall & Ulrich is good, as one would expect, except for the hep neutrinos. Since the hep and 
17 F contribute negligibly to the rate in all the detectors this is not of concern. The flux correlation 
matrix for the T c approach has all elements equal to ±1 since there is perfect correlation — i.e. all the 
fluxes depend on one parameter. This is relaxed in the updated approach including the cross section 
uncertainties [[UJ which it is claimed also reproduces table [|. 

As we shall discuss later, at least as far as flux uncertainties are concerned, the entire statistical 
content of the full solar model Monte Carlo calculation is contained in the covariance matrix V^-. 
Our method is designed to reproduce this matrix based on the matrix of flux derivatives, while the 
updated T c approach reproduces this matrix by a posteriori construction. Nevertheless, once the 
matrix is obtained from the solar models, this alone is sufficient, and there is no need for either 
approximation. For this reason, this quantity is as important to extract from solar model calculations 
as are logarithmic flux derivatives a^, and we suggest that future work on solar model calculations 
include the results for Vij explicitly. 

In our fits, to be described later, we use the updated fluxes from the Bahcall & Pinsonneault solar 



model |15j which incorporates Helium diffusion and new equation of state and opacity calculations. 
Although a full Monte Carlo treatment of the flux uncertainties has not been performed for this 
model we have updated the correlation matrix shown in table § to incorporate the errors on the 



input parameters as given in |L5[]. This does not include the uncertainty in the fluxes from variations 



in diffusion, but is the best use of currently available information. 

3. Experimental Rate Uncertainties 

The central quantity to use in determining how well model predictions agree with the observed 
rates will be the rate covariance matrix. When solar model uncertainties can be completely param- 
eterized in terms of the flux covariance matrix, the covariance matrix for the rates can be calculated 
directly from that for the fluxes. In this case, for any theoretical model the predicted rate in the 
detector R a (a =H,K,Ga) is a linear combination of the fluxes <f>j with coefficients functions of the 
theory parameters (see equation (p3|) ). If we write R a = r a j(fij with r a j = r a j(Am 2 , sin 2 26) then it is 
straightforward to show that 

Vab = ^2r aj r bk V jk (4) 

jk 

It is important to note that at this stage experimental uncertainties, including those from detection 
cross section uncertainties have not been introduced. 
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The correlation matrices for both the fluxes and the experiments (assuming standard model 
interactions for the neutrinos) are shown in tables |] and |] for the Bahcall & Ulrich standard solar 
model(s) [[nj and our power-law approach. Note that the experimental rates are almost perfectly 
correlated (a fact which was ignored in our earlier work M). The correlation between experiments 
can decrease once neutrino mixing is allowed. For example the correlations can be as low as ~ 25% 
for (Am 2 , sin 2 26) in the small-angle allowed region (see figure ^|). However generally the correlation 
is above 80%. It is initially surprising that the rates for Homestake and Kamiokande, which measure 
principally 8 B neutrinos, should be (strongly) positively correlated with Gallium, which measures 
principally pp neutrinos, when 8 B and pp neutrinos are strongly antz-correlated! The resolution of 
this apparent paradox is that while the major contribution to the Gallium rate is due to pp neutrinos, 
the 8 B and 7 Be neutrino fluxes are much more uncertain and are the principal contribution to the 
uncertainty in the Gallium rate. For Gallium 

R Ga = 71pp + 34 7 Be + 14 8 B + --- SNU (5) 

The relative errors of the pp, 7 Be and 8 B fluxes are 72%, 5% and 15% respectively. Clearly the 
uncertainty in 7 Be and 8 B dominates the uncertainty in the Gallium rate. 

4. Neutrino Transport 

In order to determine the experimental rate matrix described above, we must utilize analytic or 
numerical techniques to propagate neutrinos through the sun, empty space, and the earth in order to 
determine survival probabilities and resulting flux modulations. The methods used differ, depending 
upon whether one is interested in the region of mass and mixing angle space where MSW oscillations 
or vacuum oscillations are important. 

a. Vacuum oscillations 

In addition to the MSW model, there exists the possibility that the observed deficit of neutrinos 
could be due to oscillations in vacua between the sun and the earth |i~6| . In this case the details of 
the production in the sun are unimportant and we need keep track only of total flux variations. The 



survival probability is [12 



P{u e -^v e ) = l- sin 2 26 sin 2 — (6) 

L v 

with Ly = 4nE/Am 2 . Additionally one can average this survival probability over the change in the 
Earth-Sun distance during times comparable with the average duration of an experimental "run". 
We find our conclusions do not depend on the averaging. 

Performing a fit to the current experimental data (as described later) we find a small region in 
parameter space which is allowed at the 95% confidence level. This region agrees in general with 



those found by other authors |T8|, [17[) who have explored this theory in detail, and we will have 
nothing further to say about it. 
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b. MSW oscillations 

Perhaps the most promising neutrino mixing solution to the solar neutrino problem is the Mikheyev 
Smirnov-Wolfenstein (MSW), or matter-enhanced mixing, model |HJ. In this section we give details 
of our approximations and modelling of this effect in relation to computing the predicted rates in the 
Homestake [JB5) , Kamiokande [^0] and Gallium pl| , |22f neutrino experiments. 

To compute the rate predicted by a model for any detector we need information about the 
neutrino production in the sun. We use the flux distributions over the production regions d(f>i(r) 
and the electron number density as a function of solar radius, N e (r), from and the scale heights 
at resonance r tabulated in p3f for use with their analytic approximations^. We have explicitly 

3 /\dN e /dr 



/yres , 

e 



from the Bahcall & Pinsonneault standard solar model 15 



checked that using tq 

produces the same results. Additionally we assume that the energy spectrum of neutrinos at each r 
is as described in |[12j| . We have fitted the spectra for all species as a polynomial times the relevant 
/3-decay spectrum (correcting the typographical error in eq. 8.15 of [0]). These values are then input 
into the analytic expressions for the v e survival probability [^3[] (see also [0]), as outlined in our 
previous work |14| and summarized below. 

If the neutrino passes through a resonance on its way through the sun then we define 



Ario 




(7) 



where N e (r) is the electron density profile in the sun. The electron density at resonance is given by 

i2X / cos 29 



( Am 2 



\2E ) \V2G F J ' 



In terms of Uq we classify the transition as either adiabatic (4n ^> 1) or non-adiabatic (4n < 1). 

Let be the electron density at the point of v production, then for neutrinos in the adiabatic 
region, or those in the non-adiabatic region with iV e rcs < JV(°)/(l + tan20) the analytic expression for 
the v e survival probability is given by 



where 



1 

2 + 



X 



l + e 



1 + e- 



cos26L cos 26, 



2tt r - 



Am 2 
~2E 



y 



cos 261 



n 



2tt no(l -tan 2 6>) 

(1 - rj)/y/(l -7]) 2 + tan 2 29 



(9) 



(10) 

(11) 
(12) 
(13) 



-'^We correct a programming error in our earlier work in which ro was incorrectly read from the table. 
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For neutrinos in the non-adiabatic region produced near resonance, iVj°V(l +tan2#) < iV e res < 
/(l — tan 29), the corresponding expression for the survival probability is 



P(y e v e ) = - [1 + exp(-7rn )] , 



(14) 



while for non-adiabatic transitions with Nl es > / (l — tan 29) or adiabatic transitions with Nl es > 
we use 



P(v, 



- + - cos 29 m cos 29. 



We have also included in our analysis the effects of double resonances in the sun, see [14 



(15) 



c. Earth Effects 

It has loner been known 25 that for Am near 10" 6 eV 2 and large sin 29 it is possible to 'regen- 
erate' v e by having neutrinos pass through the earth. The survival probability P(u e — > u e ) is very 
sensitive to the path length of the neutrinos in the earth and so neutrinos with parameters in this 
range should give rise to day/night and seasonal variations in the observed flux. Since no such effect 
has been seen |20| this serves to rule out a region of parameter space near Am 2 ~ 10 _6 eV 2 and 
sin 2 29 ~ 0.2. 

We follow |26j in including this "Earth effect" in our fits, though our treatment differs from theirs. 
Rather than keep track of the predicted dependence of P{y e — * u e ) on the path length (which changes 
during the "night" ) and fit to the data in many binsQ, we choose to use a number which summarizes 
that no effect is actually seen. Consequently we use the quoted measurement of pT] 



' day — night 
day + night 



-0.08 ±0.11 



(16) 



year 



which is independent of the solar model flux uncertainties. Since this quantity does not depend on 
the neutrino flux we can simply add the x 2 from this fit to the \ 2 obtained from fitting to the time 
average rates as will be described later. The effect will be to rule out a region of parameter space 
where a large day/night effect would be predicted. 

To predict the l.h.s. of fll6l) we follow |26|, |27fl . Since only the integrated electron density along the 
line of sight matters for the average P[y e — > v e) we model the Earth as 5 concentric shells of constant 
electron density N e , which we have taken from the models of p8[ and listed in table |l|. Including the 
Earth effect the survival probability of a v e which has MSW survival probability Pmsw is given by 



27 



Pe — -PmswH 



1 — Pi 



MSW 



^2 — Pmsw 



tan29(ab* + ba* 



(17) 



2 We note in passing that the binned data of |2(| for the day/night effect has a very low \ 2 P er degree of freedom, 
which may indicate correlated (systematic) uncertainties in this data set. In any case such a low \ 2 wm bias a fit in 
which these points form most of the degrees of freedom. 
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where a and b are elements of the unitary matrix which describes the evolution of neutrinos through 
the earth 

( "e(t) )-( a b \ ( "M \ nH ) 

{ v x (t) ) {-b* a* ){ u x (0) ) [L * } 

Once we have solved for this matrix for an arbitrary shell of our model Earth we can obtain the full 
matrix by multiplication of the evolution matrices for each shell in the appropriate order (entering 
and leaving). Thus the problem reduces to calculating a and b for propagation through a shell of 
constant N e . Dropping a constant energy offset from v e and u x , which contributes only an irrelevant 
overall phase, the evolution equation is 




'G F N e Am 2 \ Am 2 



V 



where <Ji are the Pauli matrices. We solve this using the identity exp[ia- B\ = cos |o|l + isin \ a\(a ■ a 
to yield 

a = cos \h\ — ihs sin \h\ 
b = — ihisin\h\ 



(20) 



with 



'Am 2 . _ „ G F N P Am 2 



h= I— — sin 29, 0, ^1—1 — cos 26\ x path length (21) 

AE a/2 AE 




and h = h/\h\. Although our model is relatively crude, given that we are trying to fit to the absence 
of an effect it is sufficient for our purposes. 

The final task is then to integrate over the paths through the Earth during the course of the 
night/year. In our model Earth with spherical symmetry the path is totally defined by giving the 
angle #o subtended at the center of the Earth by the point of entry of the v beam and the detector. 
In the limit that the Earth-Sun distance is much larger than the Earth's radius we have that 

= sin 5 sin i + cos <5 cosi cos(0 + n) (22) 

where is the azimuthal angle between the Sun and the detector as measured from the center of 
the Earth, 5 is the detector latitude and i = 23°.5 sm(u)Qt) is the inclination of the ecliptic to the 
Earth's equator. Averaging over and u Q t we obtain the distributions for 8q, with which we can 
then determine the v e survival probability averaged over night /year. For a given mass and mixing 
angle, we compare this value with the r.h.s. of (0) in determining the x 2 fit to the data. The survival 
probability P[y e — > u e ) including the Earth effect is shown in figure [L] for a parameter set which can 
be compared with Fig. 3b of [ p7[ . 

d. Calculating the rates 
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Using the above and the formulae for Pmsw outlined in the previous section we computed the 
survival probability averaged over the night and the year. These probabilities and the fluxes for each 
species, j, are then convolved with the detector response Di{E v ) [T^, [T3|, |2T5|, |5(J for neutrinos of 
flavor i and energy E v to get the predicted rate 



R = J2 J dE » <j>j{E V )Pi{E v )Di{E v 

ij 



(23) 



for each model. We include the contributions from j = pp(10), pep(l), 7 Be(2), 8 B(30), 13 N(20) 
and 15 O(20) neutrinos, where the number in parenthesis after each species is the number of energies 
computed for each spectrum in the integration. The contribution from hep and 17 F neutrinos are less 
than 1/2% for all the experiments and can be safely ignored. Our results for the iso-SNU contours 



for the Homestake, Kamiokande and Gallium experiments compare well with those in [10 



5. Data & Model Testing 

We use the latest data for the time-averaged rate in the Homestake [IJ| , Kamiokande |2(J , SAGE 
Plfl and GALLEX [^2] experiments. Since the theory predictions for the SAGE and GALLEX 
experiments are identical and the experimental values agree within errors, we have combined the two 
rates (74 ± 20 and 79 ± 12 SNU for SAGE and GALLEX respectively) in our fit. We have added 
the statistical and systematic errors in quadrature, since they are independent. The assumption 
of a gaussian distribution for the systematic error is problematic because, by its very nature, the 
systematic error has no statistical distribution. However, if we regard the gaussian as representing 
the state of our knowledge about the systematic error then clearly a function which penalizes large 
'errors' is more appropriate than a flat distribution (which would correspond to maximal ignorance 
of the size of the systematic error). We have chosen a gaussian for simplicity. The three experimental 
rates, as a fraction of the standard solar model predictions, are shown in table 2. In our analysis we 



have added the cross section uncertainties JTT| in quadrature to the quoted errors. For Gallium this 
ignores the energy dependence of the uncertainty from the resonance, but this uncertainty affects 
primarily the 8 B contribution to the rate which is already small and which we expect to be suppressed 
for the masses and mixing angles of interest to us. 

We have used two parametric methods: a \ 2 goodness-of-fit procedure and a Bayesian likelihood 
analysis [|13|] to compare the measured rates R™±<r a (a=H,K,Ga) to the rates predicted by the model 
(Am 2 , sin 2 29). Both methods rely on the assumption that the errors in the solar model predictions 
are gaussian under small variations in the solar model input parameters, which appears to be a good 
assumption |12[. For more discussion of the methods we use see [O, pTl EO, B3]. 



Including the solar model uncertainties, each set of MSW parameters (Am 2 , sin 2 29) defines a 
distribution of rate triplets R a . One can calculate the covariance matrix V^ b M for the triplets analo- 
gously to equations (@,|). To the theoretical solar model covariance, Vf b M , we add the experimental 
errors, a a , to obtain the full covariance matrix: V c 



ab 



ySM , 2 s 
v ab + <J n O, 



ab ■ 



Defining [13| 



X 



(R a — R a )V ab l (R h — R b ), 



(24) 
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where R is the standard solar + neutrino-mixing model rate prediction, the distribution of x 2 s defined 
by the theory is a chi-square distribution with 3 degrees of freedom. The statement that the theory 
is ruled out at some confidence level is the claim that the x 2 f° r the measured triplet lies in the 
large-% 2 tail of the distribution defined by the theory, i.e. that the measured value is unlikely. If the 
data is within the 95% confidence level of a given theory we say the theory parameters are allowed 
at the 95% confidence level by the data. 

Previous authors have generally implemented instead a best-fit procedure that attempts to esti- 
mate the parameters Am 2 and sin 2 29 from the data and assign errors to the inferred values. One 
takes the parameters which minimize \ 2 as the central values, with an allowed range given by the 
condition that x 2 (Am 2 , sin 2 2$) < xLin + v i with v determined by the range of a desired and the 
number of parameters being estimated. Such an approach is based on the maximum likelihood pro- 
cedure under the assumption that the correlation matrix is independent of the parameters being 
estimated. (This assumption is obviously not true for this case, but the errors introduced turn out to 
be numerically small.) This approach makes the additional assumption that x 2 (Am 2 , sin 2 29) is well 
approximated by a quadratic over the relevant range of parameters. As can be seen in figure |2| this 
assumption is clearly false over the range of (Am 2 , sin 2 29) of interest. It is important to realize that 
the statistical answers one gets depend upon the questions one asks! This method does not address 
the question of what regions of model space are allowed by the data, but rather what regions provide 
a best fit under the assumption that the model is correct, for some set of parameters. The allowed 
region determined differs from that for the method outlined above as it asks a different statistical 
question: not what models are allowed by the data but what are the errors on the best-fit Am 2 and 
sin 2 29. 

In addition, an approach for calculating allowed regions has recently been advocated which 
uses non-standard definition of x 2 ■ m comparing their method to solar model Monte Carlos [11] the 
authors define "x 2 " in terms of the logarithm of the "average probability" rather than computing V a \, 
for the Bahcall & Ulrich solar models directly and using equation fl24j) . The distribution of this "x 2 " 
will not be chi-squared, and will not take into account correlations in the rates in a well defined way. 
To consistently use such a statistic, the correct distribution and confidence levels to be associated 
with it would need to be calculated. 

An alternative method, which is similar in spirit to the best fit approach, is to calculate the 
2D likelihood function £(Am 2 , sin 2 29), again under the assumption that the variations in model 
predictions (and the experimental errors) are gaussian. In this case the likelihood function is defined 



C oc , exp 
VdetV 



I 2 

2 X 



(25) 



3 Note that the authors of [[10| define a likelihood function as a sum of gaussians and redefine v above to give 
regions consistent with this approximate likelihood function. While the statistical meaning of this hybrid method is 
not immediately clear, the final regions obtained are not much different than provided by more conventional statistical 
treatments. 
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By use of Bayes' theorem, the conditional probability for Am 2 and sin 2 26, given the experimental 
measurements, is proportional to the likelihood function times the a priori probability distribution 
for the MSW parameters (for an alternative interpretation see p3|), which is usually referred to 



as the posterior distribution. If we assume, from scaling arguments, that logarithmic intervals in Am 2 
and sin 2 26 are equally likely, before any experiment is performed, then the posterior distribution is 
simply proportional to the likelihood function £(log Am 2 , log sin 2 26). 

To calculate the 95% confidence regions for Am 2 and sin 2 26 we follow the method used in assigning 
regions for gaussian distributions (which C is not). Let us define a region in parameter space 

A(X) = { (log Am 2 , log sin 2 26) | £(log Am 2 , log sin 2 26) > Xj . (26) 

Then 

r(A)= / £(logAm 2 ,logsin 2 2#) d(log Am 2 )d(logsin 2 20) (27) 

■>A(X) 

is a continuous, monotonic decreasing function of A, and the 95% confidence region is given by 
A(X*), where r(A*) = 0.95r(0). (We note that this method is somewhat arbitrary for multiply- 
peaked likelihood functions such as ours, but it is nonetheless well defined.) This confidence region 
is interpreted as the region that contains, with 95% probability, the true values of Am 2 and sin 2 26. 
Although the interpretation of the region is different than that allowed by the x 2 method, the two 
regions are encouragingly similar. In the limit that the likelihood function were gaussian (x 2 is a 
quadratic function of Am 2 and sin 2 26 and det V is constant) the regions would be ellipses as in ||10|| . 
Thus the departure from elliptical shape is an indication that the likelihood function is not simply 
gaussian. 

6. Results 

Our principal result, the 95% C.L. allowed regions in MSW parameter space, based only on 
statistical uncertainties in the present formulation of the standard solar model, is shown in figure 
H] for both the x 2 an d £ methods. The regions shown are obtained by requiring that x 2 < 9-49 (4 
dof), including both the rate and the day/night fits. We also show the region obtained by requiring 
X 2 < Xmin + 6.0 (2 parameters) for comparison. 

We see that there are two allowed regions, a large mixing angle (adiabatic) region and a small 
mixing angle (non-adiabatic) region. The small mixing angle region is favoured over the large mixing 
angle region, though both are "allowed" at the 95% confidence level. Using the likelihood function 
we can ask what are the relative probabilities of the large and small angle regions, e.g. we find 
P(sin 2 26 > 0.1) ~ 0.3P(sin 2 26 < 0.1) (see |I| for a different way to ask this question). 



One of the largest uncertainties in calculating the expected rates comes from S\j, the nuclear 
cross section parameter for the reaction 7 Be(p, 7) 8 B. This uncertainty directly affects the flux of 
8 B neutrinos, and is due to both experimental uncertainties and the difficulty of extrapolating the 
experimental results to the low energies relevant in the solar interior. There is a significant difference 
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between the cross sections inferred from the two experiments which have been performed at the lowest 

use a 



energies, indicating some significant systematic error in this parameter. The authors of ref |L5 
value for Sij that is intermediate to these two results, with errors which do not overlap the central 
values. This value and its errors, which we have used in determining the correlation matrix, may 
not properly reflect the uncertainty in Sij. In order to explore the implications of a larger estimate 
for the error in S\j, in figure ^ we present the allowed regions assuming as the error for S\j the 
difference between the two central values of ref 35], which corresponds to a 21% uncertainty. As 



is expected the allowed regions are correspondingly increased. 

Recent refinements in solar models, i.e. including heavy element diffusion, have changed the 
predicted neutrinos fluxes from those of the Bahcall & Pinsonneault model. As we have indicated, 
changes such as this, which includes new physics rather than new numerical values for the input 
parameters and their errors can have a large effect on the allowed regions. The situation with respect 
to these new solar models is still not settled, and the correlation matrix including parameters for 
the effects of heavy element and helium diffusion are not yet available. Nevertheless, to estimate the 
magnitude of the effect of such changes we have use a hybrid procedure which uses neutrino fluxes 
including heavy element and helium diffusion from [37. 38j, but our old correlation matrix. Specifically 



we have artificially changed the fluxes of the Bahcall & Pinsonneault model by percentages equal to 
those shown in table |5] but used the flux correlation matrix from table |j. While this method is not 
fully consistent, and the fluxes used are preliminary, it should approximate the main effects of these 
changes and illustrate the possible shift in the allowed regions that can be expected for such models. 

We display in figure f| our result for the allowed range of parameter space in this case. The change 
in the allowed region reinforces our earlier remarks: the potential change in the allowed regions due 
to such modifications of the solar model can be larger than indicated by the inclusion of the usual 
solar model uncertainties. Thus it is prudent to realize that the presently "allowed" regions are now 
only suggestive. Further solar model improvements could change their shape, and position. 

7. MSW and Refined Confidence Limits: Future Work 

The statistical analysis we have provided here is straightforward, and resolves various inconsis- 
tencies present in previous analyses. As such, it should provide a firm basis with which to analyse 
future results. However, for the MSW solution, neither this approach, nor any other to date actually 
properly accounts for all solar model uncertainties when determining allowed ranges in parameter 
space. As we have just indicated, statistical solar model uncertainties do not incorporate possible 
systematic shifts in fluxes due to new, non-exotic, physics, which could dramatically alter the shape 
of allowed regions. Beyond this, however, in the case of the MSW solution no set of neutrino flux 
uncertainties can carry all of the relevant information on solar model variations. This is because 
the neutrino deficits which result from traversing the solar interior themselves depend sensitively 
on the details of the solar density and temperature. Thus, simply calculating the initial neutrino 
fluxes over a wide range of solar models, but propagating them using the density-temperature re- 
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lationship of the standard solar model is not fully consistent. Stated in a language which can be 
compared with that given above, when determining predicted rates in detectors R a , we can no longer 
write R a = r a j(pj with r a j = r a j(Am 2 , sin 2 29). Rather, r a j now becomes also a function of N e (r), 
r a j = r nj (Am 2 , sin 2 29, N e ), so that the decomposition which lead to equation (Q) can no longer be 
carried out. Instead, for each value of Am 2 and sin 2 29 one must carry out the average over solar 
models directly, accounting for propagation in the sun in each model, in order to determine the 
proper predicted rate covariance matrix. We do not know a priori how large an effect such a detailed 
accounting will produce, although there is some reason to expect it will not change the allowed region 
drastically. Computationally this is far more daunting, and we will report on this analysis in a future 
publication |3(J where we will carry out such a procedure for the newest set of solar models of Bahcall 
& Pinnsonault, in which heavy metal diffusion is accounted for, and for which the predicted neutrino 
fluxes appear to be somewhat larger. 

8. Future Experiments 

At very large mixing angles we expect that the v e survival probability will be roughly independent 
of energy so that the spectrum of neutrinos seen would be unchanged but for the normalization. 
In the adiabatic region the existence of a resonance implies a large suppression of the v e survival 



probability while the converse is true in the non-adiabatic region ||12[ . Hence we expect that for 
(Am 2 , sin 2 29) in the small mixing angle allowed region, the lower energy pp and 7 Be neutrinos, 
which have energies that correspond to the adiabatic regime, will be preferentially depleted, while 
the higher energy 8 B neutrinos have a higher survival probability due to nonadiabatic level jumping. 
Two neutrino experiments currently under construction, the Sudbury Neutrino Observatory (SNO) 
and SuperKamiokande, may have the ability to detect the distortion in the 8 B neutrino energy 
spectrum for small angle MSW solutions. For most of the small angle region, SNO should be able to 
discern the spectral distortion, while it is very unlikely that the minimal shape distortion produced 
by MSW parameters in the large angle region could be detected. 

SNO of course can also measure the ratio of charged current events to neutral current events, which 
provides an indicator for v e oscillations into another active neutrino species. A ratio significantly less 
than that expected for the SSM (i.e. the ratio of the electron neutrino charged current to neutral 
currents cross sections [0) would be a strong indication of MSW mixing. However, the neutral 



current events are signaled by the production of a free neutron, and the background for this process 
can be problematic. 

The improved statistics of SuperKamiokande, which expects to see on the order of 8000 events per 
year JIT], can be used to examine more closely the effects of v e regeneration through the earth (see 
section 4c). In figure [5|, we plot the contours for several values of (day-night) /(day+night), superim- 
posed upon the allowed regions. Recall that our model for calculating the above ratio automatically 
assumes an average of the year and the entire night and uses a very crude model of the earth. An 
analysis using the methods of section 4c on a more realistic model of the earth could be performed, 
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however figure [5] serves to show the potential for narrowing the allowed regions if a positive day/night 
variation is seen. If SuperKamiokande does find a signal for day/night variation then binning the 
data vs cos<5 SU n (c.f. eq. ^2] with 5 sun = (9 + tt)/2) would provide more information than our simple 
average. 

An experiment sensitive to 7 Be neutrinos can potentially discriminate between the large and small 
mixing angle solutions in addition to confirming the depletion of the 7 Be flux. The predicted rate for 
Borexino [42| as a fraction of the standard solar model rate is [43] 



^Borexino = 0.787P(z/ e -> u e ; 7 Be) + 0.213 (28) 

where the 0.213 is absent for oscillations into a sterile neutrino. For the two currently allowed MSW 
parameter regions, the predicted rates in Borexino are shown in figure A Borexino rate of less than 
0.3 would provide not only a striking confirmation of the solar neutrino problem, but also indicate 
the small mixing angle region of the MSW solution, whereas rates between 0.5 and 0.8 would point to 
the large mixing angle solutions. A detected rate of about 0.35 of the standard solar model would not 
allow discrimination between the two solutions, but would nonetheless be further evidence in support 
of new neutrino physics. However note that after lyr of running Borexino can at best measure the 
rate to ~ 30%. 

9. Implications for particle physics models 

If future experiments confirm the deficit of electron neutrinos indicated by the current data, we 
would have the first (indirect) evidence of physics beyond the standard model: neutrino masses. 
Further solar neutrino studies coupled with upcoming neutrino oscillation experiments [H are the 



current best hope of seeing neutrino masses in the cosmologically interesting range J2 m i — 3 — 30eV. 
(The region of mass-mixing angle space of interest for oscillations which may explain the deficit of 
atmospheric muon neutrinos can be probed by several proposed long baseline oscillation experiments 
f45fl.) The mass- mixing angle parameters implied by the allowed regions shown in figure |2|, while not 



indicative of any particular particle physics models for neutrino masses, are consistent with models 
which incorporate a seesaw mechanism. Many of these models can also accommodate the observed 



deficit of the ratio of atmospheric v^j v e |46| and in some cases also allow for the v T to be cosmological 



hot dark matter |47| or provide contributions to neutrinoless double /3-decay at the level of m ~ leV 



There is a significant literature in the particle physics community on constructing models which go 
beyond the Standard Model of Electroweak interactions and many of these models have interesting 
implications for neutrino properties. Here we discuss some classes of models which are currently 
popular and which relate directly to the solar neutrino problem. 

Models in which the neutrino mixing angles are similar to the CKM angles in the quark sector 
f49|| now appear to be disfavoured by the data. A class of models based on grand unification particle 



physics models and a see-saw mechanism for neutrino masses ED[ give masses and mixings which 
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can lie in the small angle region of figure |5"1| , |52| . In some cases these models can also incorporate 
solutions to the atmospheric neutrino deficit, provide the hot dark matter component of currently 
popular mixed dark matter models (on the order of 1-10 eV in neutrino mass), and even accommodate 



a Majorana mass of 1-2 eV for neutrinoless double (3 decay [p3 |. In these models the masses of the 
light neutrinos we see are a combination of the Dirac masses of the usual neutrinos plus new right 
handed neutrinos vr which additionally have Majorana masses. The ur are placed in GUT gauge 
group multiplets along with the quarks and leptons and get masses from the same Higgses. This 
relates the Dirac mass matrices of the neutrinos in these models to those of the quarks and leptons. If 
further "textures" are assumed for the heavy Majorana mass matrix of the z/r, one obtains predictions 
for the masses and mixings of the observed light neutrinos, usually with one free (overall mass) scale 
and a small number of group theory factors. We note in passing that the presence of (powers of) 
these group theory factors can significantly alter the naive see-saw predictions. 



Other models exist |5J, [55], |56|, [57J which generate masses and mixing angles in the large angle 
allowed region of ^. Such models can allow for simultaneous solution of the solar and atmospheric 
neutrinos [5(| or link solar neutrinos with double /3-decay experiments |55[ or both |58| . 



Some authors have considered a radiative mechanism for the generation of neutrino masses, and 
found models which can accomodate two of the three neutrino mass solutions (solar neutrinos, at- 
mospheric neutrino deficit, dark matter) [ |6C|j . 

Further constraints on models for neutrino masses which invoke oscillations into sterile neutrinos v s 
are obtained by considering big bang nucleosynthesis |JT|. The large angle region is excluded for v e —v s 
oscillations based on present observations of the primordial 4 He abundance. This also eliminates 

— v s solutions to the atmospheric deficit. Arguments derived from supernova considerations 
can also be used to constrain oscillations [62], 63| . For sterile neutrinos the region of mass- mixing angle 



space restricted by these arguments, while of interest for sterile neutrinos as dark matter candidates, 
is not relevant for solar neutrino oscillations [63 . 



10. Conclusions 

In this paper we have presented an updated analysis of the implications of the four currently 
operating solar neutrino experiments. Our analysis incorporates a straightforward and comprehensive 
treatment of the known theoretical statistical uncertainties, which we have outlined in detail. We 
have given a full account of our methods and assumptions so others can compare with our work. 

We find that the current solar neutrino experiments provide a useful constraint on the masses 
and mixing angles of neutrino in models where neutrino mixing is the resolution of the solar neutrino 
problem. All the quantifiable errors in established solar models are included in this constraint. Both 
resonant (MSW) and non-resonant (just-so) neutrino oscillation models are allowed by the data. In 
considering the implications of these figures, it is important to note that systematic uncertainties 
remaining in both the solar model calculations and input parameters can have an effect on the 
properties of the allowed regions, as shown in figure 0L Nevertheless, as solar models improve, the 
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consistent statistical analysis we have denned here will continue to gain in significance. 

Future experiments have great promise for confirming the solar neutrino problem and firmly 
establishing the need for neutrino based solutions. In particular we have examined the potential of 
SuperKamiokande, SNO and Borexino to provide further constraints on the masses and mixing angles 
of neutrinos in such models. If the charged to neutral current ratio measured by SNO indicates the 
probability of neutrino oscillations, the presence (absence) of spectral distortion will further constrain 
the mixing parameters to the small (large) angle regions. (While SNO is not capable of distinguishing 
between large angle oscillations into sterile neutrinos and solar model solutions, such oscillations have 
already been ruled out by big bang nucleosynthesis as mentioned above |61].) Borexino also possesses 



the ability to distinguish between small and large angle MSW regions to some extent. If the large 
angle solution turns out to be favoured, then SuperKamiokande should provide a sensitive probe of 
the allowed mass and mixing through the measurement of the day/night effect. 
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R Rca 


p 


(V2G F N e ) 


0.0000-0.1910 


12.858 


4.87 


0.1910-0.5471 


11.024 


4.18 


0.5471-0.8948 


4.964 


1.88 


0.8948-0.9341 


3.923 


1.49 


0.9341-1.0000 


2.292 


0.87 



Table 1: The model Earth that we used in calculating the regeneration effect. Densities in column 2 
are given in g/cm 3 and the final column is in 10~ 7 eV 2 /MeV, assuming n p = n n for the Earth interior. 



Experiment 


Rate 


Homestake 

Kamiokande 

Gallium 


0.31 ±0.03 
0.51 ±0.07 
0.59 ±0.08 



Table 2: The experimental rates, normalized to the standard solar model predictions, used in the fits. 
The rates for SAGE and GALLEX have been combined and cross section uncertainties for Homestake 
and Gallium have been added, in quadrature, to the experimental errors. 
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Table 3: The experiment and flux correlations (xlOO) computed using the 1,000 solar models of 
Bahcall & Ulrich. 
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Table 4: The experiment and flux correlations (xlOO) computed using the power-law Monte Carlo 
approach. 



Flux 


% Change 


pp 


n 


pep 


n 


hep 


n 


7 Be 


5| 


8 B 


14T 


13 N 


4| 


15Q 


24T 



Table 5: The percentage change in the fluxes for each neutrino species in the Profntt solar model 
arising from including heavy element diffusion. 



19 



Figure 1: The survival probability P{y e — > v e ) including the Earth effect for mixing angle sin 29 = 0.4. 
The dotted line is the MSW probability without the Earth effect, the solid line is the probability for 
a neutrino passing through the center of the Earth and the dot-dashed line shows the probability 
averaged over the night/year. 



Figure 2: Region of MSW mass-mixing angle space allowed at the 95% confidence level for the com- 
bined Homestake-Kamiokande-Gallium data including solar model uncertainties from the x 2 analysis 
(solid). Also plotted is the 95% confidence region from the likelihood function analysis (dotted) and 
the region obtained by requiring x 2 < Xmin + v (dashed). 



Figure 3: Likelihood function £(log Am 2 , log sin 2 29) for the combined Homestake-Kamiokande- 
Gallium data. 
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Figure 4: Region of MSW mass-mixing angle space allowed at the 95% confidence level in the Bahcall 
& Pinsonneault standard model for the combined Homestake-Kamiokande-Gallium data including 
normal solar model uncertainties for x 2 < 9.49 (solid). Also plotted is the 95% confidence region 
increasing the error on S±j to 21% (dotted) and the region obtained for fluxes approximating the 
effects of metal diffusion (dashed, see text). 



Figure 5: The allowed region from the fit to the experimental rates, plus the contours of (day- 
night)/(day+night) rate. The contours are 1% (solid), 5% (dotted), 10% (dashed) and 15% (long- 
dashed). An average over the night and the year has been assumed for these contours, and the model 
of the interior of the earth used was very simplistic. 



Figure 6: Predictions for the Borexino event rate as a fraction of the standard solar model value. 
The histograms show relative frequencies of predicted event rates in the large (dashed) and small 
(solid) angle regions. 
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